﻿#include <stdio.h>
#include <string.h>
#include <math.h>
#include<stdlib.h>

#define M 64
#define T 16
int    Prec=2;        // Setting the precision of the answer.
double Appr=1E-6;  // Using a very tiny number instead of lim{x->0}.
double Inf=1E8;    // Using a very large number instead of infinity.
int    Inf_r=1024; // Using a very large number instead of infinity in series.
// Using operators instead of functions.
//Expressions
#define sin_e  "1s"
#define cos_e  "1c"
#define tan_e  "1t"
#define asin_e "1S"
#define acos_e "1C"
#define atan_e "1T"
#define sinh_e "1i"
#define cosh_e "1o"
#define tanh_e "1a"
#define exp_e "1p"
#define ln_e  "1n"
#define lg_e   "1g"

//Operators
#define sin_o  's'
#define cos_o  'c'
#define tan_o  't'
#define asin_o 's'
#define acos_o 'C'
#define atan_o 'T'
#define sinh_o 'i'
#define cosh_o 'o'
#define tanh_o 'a'
#define exp_o  'p'
#define ln_o   'n'
#define lg_o   'g'

typedef enum {
    normal,oper_mis, brac_mis, beyond_range, lim_diffuse
} Err;
Err cur_err=normal;
// stands for operators mismatching,brackets mismatching,out of definition.

//Interface
double cal(const char *equ_argu);
double fx(const char *equ_argu, char *var_str,double *val_arr);
double fx1(const char *equ_argu, double x);
double diff(const char *equ_argu, char *var_str, double *val_arr, int n);// diff for var[0]
double integral(const char *equ_argu,double lower, double upper);
double recur(const char *equ_argu, double *val, int val_i, int n);
double sumrec(const char *equ_argu, double *val, int val_i, int n);
double limrec(const char *equ_argu, double *val, int val_i);

//Test
void test_cal(void) {
    char equ[] = "sin(1+1)+1";
    printf("%s = %.2f\n",equ,cal(equ));
}

void test_func(void) {
    char equ[] = "sin(x*y)+1",equ2[] = "sinx";
    char var[] = "xy";
    double x = 2,val[] = {1,2};
    printf("\nf(x,y) = %s\n",equ);
    printf("f(%.2f,%.2f) = %.2f\n",val[0],val[1],fx(equ,var,val));
    printf("f'x(%.2f,%.2f) = %.2f\n",val[0],val[1],diff(equ,var,val,1));
    printf("f''x(%.2f,%.2f) = %.2f\n",val[0],val[1],diff(equ,var,val,2));
    printf("\ng(x) = %s\n",equ2);
    printf("g(%.2f) = %.2f\n",x,fx1(equ2,x));
    printf("{%.2f,%.2f}int g(x) = %.2f\n",x,x+1,integral(equ2,x,x+1));
}

int test_rec(void) {
    char equ[] = "a[n-1] + a[n-2]";
    double val[] = {1,1},val_i = 2;
    int n = 100;
    printf("\na[n] = %s\n",equ);
    printf("a[%d] = %.2f\n",n,recur(equ,val,val_i,n));
    printf("s[%d] = %.2f\n",n,sumrec(equ,val,val_i,n));
}

int main(void)
{
    test_cal();
    test_func();
    test_rec();
    return 0;
}

inline int isoper(char c) {
    char opers[]={'^','*','/','+','-',sin_o,cos_o,tan_o,asin_o,acos_o,atan_o
    ,sinh_o,cosh_o,tanh_o,exp_o,ln_o,lg_o,'\0'};
    if(strchr(opers,c) )
         return 1;
    else
        return 0;
}

char *replace(char *beg, char *end,char *src) { //replace [beg,end) with string src.
     char temp[M];
     strcpy(temp,src);
     strcpy(temp+strlen(src),end);
     strcpy(beg,temp);
     return beg;
}

char *remove_ch(char *ch) {
    int i;
    for (i=0;ch[i];i++)
        ch[i]=ch[i+1];
    return ch;
}

double *remove_val(double *arr,int index, int len) {
     int i,l;
     l=len-index;
     for (i=index;i<=l;i++) {
        arr[i]=arr[i+1];
     }
     return arr;
}

char *pretreat(char *equ) {
   //Turn f(x) into equation like "a@x",
   //in which "a" is a formal lvalue,"@" is a operator in place of function,x is the variable,
   //so that cal_poly can treat f(x) as other equations such as "1+2".
   int i;
   for (i=0;equ[i];i++) {
     if (equ[i]== 's' && equ[i+1]=='i' && equ[i+2] == 'n')
        replace(equ+i,equ+i+3,sin_e);
    else if (equ[i]== 'c' && equ[i+1]=='o' && equ[i+2] == 's')
         replace(equ+i,equ+i+3,cos_e);
    else if (equ[i]== 't' && equ[i+1]=='a' && equ[i+2] == 'n')
        replace(equ+i,equ+i+3,tan_e);
    else if (equ[i]== 'a' && equ[i+1]=='s' && equ[i+2] == 'i' && equ[i+3] == 'n')
        replace(equ+i,equ+i+3,asin_e);
    else if (equ[i]== 'a' && equ[i+1]=='c' && equ[i+2] == 'o' && equ[i+3] == 's')
        replace(equ+i,equ+i+3,acos_e);
    else if (equ[i]== 'a' && equ[i+1]=='t' && equ[i+2] == 'a' && equ[i+3] == 'n')
        replace(equ+i,equ+i+3,atan_e);
    else if (equ[i]== 's' && equ[i+1]=='i' && equ[i+2] == 'n' && equ[i+3] == 'h')
        replace(equ+i,equ+i+3,sinh_e);
    else if (equ[i]== 'c' && equ[i+1]=='o' && equ[i+2] == 's' && equ[i+3] == 'h')
        replace(equ+i,equ+i+3,cosh_e);
    else if (equ[i]== 't' && equ[i+1]=='a' && equ[i+2] == 'n' && equ[i+3] == 'h')
        replace(equ+i,equ+i+3,tanh_e);
    else if (equ[i]== 'e' && equ[i+1]=='x' && equ[i+2] == 'p')
        replace(equ+i,equ+i+3,exp_e);
    else if (equ[i]== 'l' && equ[i+1]=='g')
        replace(equ+i,equ+i+3,lg_e);
    else if (equ[i]== 'l' && equ[i+1]=='n')
        replace(equ+i,equ+i+3,ln_e);
   }
   return equ;
}

char *subs(char *equ ,char *var_str,double *val_arr) {
    int equ_i,var_i=0;
    char temp[M]={'\0'};
    pretreat(equ);
    for (var_i=0;var_str[var_i];var_i++) {
      sprintf(temp,"%lf",val_arr[var_i]);
      for (equ_i=0;equ[equ_i];equ_i++) {
         if (equ[equ_i] == var_str[var_i]) {
            replace(equ+equ_i,equ+equ_i+1,temp);
         }
      }
    }
    return equ;
}

double sum_t (double *val_arr, int l) {
   int i;
   double sum=0;
   for (i=1;i<=l;i++) {
       sum += val_arr[i];
   }
   return sum;
}

double cal_poly(const char *equ_argu) {
    char oper_que[M]={'\0'}; // a string consists of operator
    char temp[M] = {'\0'};
    double val_que[M]={0};   // a queue consists of value
    int beg=0,end,oper_i=0,val_i=0;
    char equ[M];
    strcpy(equ,equ_argu);
    pretreat(equ);
    for (end=1;equ[end];end++) {
        if (isoper(equ[end])) {
            //push the val into queue
            strncpy(temp,equ+beg,end-beg);
            temp[end-beg]='\0';
            sscanf(temp,"%lf",&val_que[val_i]);
            val_i++;
            //push the oper into the queue
            oper_que[oper_i]=equ[end];
            oper_i++;
            //reset beg,end
            beg=end+1;
            end++;
            //end++,for(;;end++) => end += 2,
            //If there is a negative sign, it will regard it as part of the value
            //without judging whether it is an operator.
        }
    }
     strncpy(temp,equ+beg,end-beg);
     temp[end-beg]='\0';
     sscanf(temp,"%lf",&val_que[val_i]);
     val_i++;       // val_i represents the len of val_arr.
     oper_que[oper_i]='\0';
     // Calculate an expression and replace it with its answer.
     // The operator, oper_que[i], has two operands, val_que[i] and val_que[i+1].
     //priority
      for (beg=0;oper_que[beg];beg++) {
        if (oper_que[beg] == sin_o) {
            val_que[beg+1]=sin(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
        else if (oper_que[beg] == cos_o) {
            val_que[beg+1]=cos(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == tan_o) {
            val_que[beg+1]=tan(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == asin_o) {
            val_que[beg+1]=asin(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == acos_o) {
            val_que[beg+1]=acos(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == atan_o) {
            val_que[beg+1]=atan(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == sinh_o) {
            val_que[beg+1]=sinh(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == cosh_o) {
            val_que[beg+1]=cosh(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == tanh_o) {
            val_que[beg+1]=tanh(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == exp_o) {
            val_que[beg+1]=exp(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == ln_o) {
            val_que[beg+1]=log(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
         else if (oper_que[beg] == lg_o) {
            val_que[beg+1]=log10(val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
     }
     for (beg=0;oper_que[beg];beg++) {
        if (oper_que[beg] == '^') {
            val_que[beg+1]=pow(val_que[beg],val_que[beg+1]);
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
     }
     for (beg=0;oper_que[beg];beg++) {
        if (oper_que[beg] == '*') {
            val_que[beg+1]=val_que[beg]*val_que[beg+1];
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
        else if (oper_que[beg] == '/') {
            val_que[beg+1]=val_que[beg]/val_que[beg+1];
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
     }
     for (beg=0;oper_que[beg];beg++) {
        if (oper_que[beg] == '+') {
            val_que[beg+1]=val_que[beg]+val_que[beg+1];
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
        else if (oper_que[beg] == '-') {
            val_que[beg+1]=val_que[beg]-val_que[beg+1];
            remove_val(val_que,beg,val_i);
            remove_ch(oper_que+beg);
            val_i--;
        }
     }
     return val_que[0];
}

double cal_brac(const char *equ_argu) {
   int brac_stack[M]={0};
   int brac_i=0,equ_i;  //brac_i stands for the next index avaliable.
   double ans;
   char temp[M]={'\0'};
   char equ[M];
   strcpy(equ,equ_argu);
   pretreat(equ);
   for (equ_i=0;equ[equ_i];equ_i++) {
       if (equ[equ_i] == '(') {
            brac_stack[brac_i]=equ_i;
            brac_i++;
        }
        else if (equ[equ_i] == ')') {                   //error transmit wrong equation
            equ[equ_i]='\0';
            ans=cal_poly(equ+brac_stack[brac_i-1]+1);
            sprintf(temp,"%lf",ans);
            replace(equ+brac_stack[brac_i-1],equ+equ_i+1,temp);
            brac_i--;
        }
   }
   ans=cal_poly(equ);
   return ans;
}

double cal(const char *equ_argu) {
    return cal_brac(equ_argu);
}

double fx(const char *equ_argu, char *var_str,double *val_arr) {  //function of several variables
     char equ[M];
     strcpy(equ,equ_argu);
     subs(equ,var_str,val_arr);
     return cal(equ);
}

double fx1(const char *equ_argu, double x)  { //function of one variable
    char equ[M];
    double val_arr[1];
    val_arr[0]=x;
    strcpy(equ,equ_argu);
    subs(equ,"x",val_arr);
    return cal(equ);
}

double diff_sec(const char *equ_argu, char *var_str, double *val_arr, int n) { //
  double y1, y2, dy,val_temp[M];
  int i;
  char equ[M];
  for (i=0;var_str[i];i++) {
    val_temp[i]=val_arr[i];
  }
  strcpy(equ,equ_argu);
  if (n == 1) {
    y1=fx(equ,var_str,val_temp);
    val_temp[0] -= Appr;      // val_temp => x+dx
    y2=fx(equ,var_str,val_temp);
    dy = (y1 - y2) / Appr;
  }
  else {
    dy = (diff_sec(equ,var_str,val_temp, n - 1) - diff_sec(equ,var_str,val_arr,n-1)) / Appr;
  }
  return dy;
}

double limfx (const char *equ_argu, char *var_str,double *val_arr ) {
    //If f(a)-f(b)= o(a-b), we regard f(a) as the limit.
     int i;
     double fa,fb,d,val_temp[M];
     char equ[M];
     strcpy(equ,equ_argu);
     for (i=0;var_str[i];i++) {
        val_temp[i]=val_arr[i];
     }
     fa= fx(equ,var_str,val_arr);
     val_temp[0] += Appr;
     fb= fx(equ,var_str,val_temp);
     d=(fa>fb) ? (fa-fb) : (fb-fa);
     if (d <= Appr) {
        return fa;
     }
     else {
        cur_err=lim_diffuse;
        return -1;
     }
}

double romberg(const char *equ, double lower, double upper)
{
    int i, k, m, n;
    double h, x, s, p, q, ep;
    double y[T];
    //initial
    h = upper - lower;
    y[0] = h*(fx1(equ,lower) + fx1(equ,upper))/2.0;
    m = 1;
    n = 1;
    ep = Appr + 1.0;
    //iteration
    while ((ep >= Appr) && (m < T))
    {
        p = 0.0;
        for (i=0; i<n; i++)
        {
            x = lower + (i+0.5)*h;
            p = p + fx1(equ,x);
        }
        p = (y[0] + h*p)/2.0;
        s = 1.0;
        for (k=1; k<=m; k++)
        {
            s = 4.0*s;
            q = (s*p - y[k-1])/(s - 1.0);
            y[k-1] = p;
            p = q;
        }

        p = fabs(q - y[m-1]);
        m = m + 1;
        y[m-1] = q;
        n = n + n; h = h/2.0;
    }

    return q;
}

double diff(const char *equ_argu, char *var_str, double *val_arr, int n) {
       return diff_sec(equ_argu,var_str,val_arr,n);
}
double integral(const char *equ_argu,double lower, double upper) {
    char equ[M];
    strcpy(equ,equ_argu);
    return romberg(equ,lower,upper);
}

double recur_1(const char *equ, double x1, int n) {
  int i;
  for (i = 0; i < n - 1; i++) {
    x1 = fx1(equ,x1);
  }
  return x1;
};

double limrec_1(const char *equ, double x1) {
   double y,z,d;
   y=recur_1(equ,x1,Inf_r);
   z=recur_1(equ,x1,Inf_r-1);
   d=(y>z)?(y-z):(z-y);
   if (d < (double)1000*Appr) {
    return y;
   }
   else {
     cur_err=lim_diffuse;
     return -1;
   }
}

double sumrec_1(const char *equ, double x, int n) {
  int i;
  double s = 0;
  for (i = 0; i <n; i++) {
    x = recur_1(equ, x, i+1);
    s += x;
  }
  return s;
}

double recur_advan(const char *equ_argu, double *val_arr, int val_i, int n) {
  // val_i out of the end,
  // val_arr begin at val_arr[1] and is long engouh for n.
  char equ[M], temp[M];
  int equ_i, brac_i, off, i;
  double sum = 0;
  for (val_i; val_i <= n; val_i++) {
    strcpy(equ, equ_argu);
    for (equ_i = 0; equ[equ_i]; equ_i++) {
      if (equ[equ_i] == 'a' || equ[equ_i] == 's') {
        // get off_set
        for (brac_i = equ_i; equ[brac_i] != ']' && equ[brac_i]; brac_i++);
        if (!equ[brac_i]) {
          cur_err=brac_mis;
          return -1;
        }
        strncpy(temp, equ + equ_i + 4, brac_i - equ_i - 4);
        sscanf(temp, "%d", &off);
        // get index val_i-off
        if (equ[equ_i] == 'a') {
          sprintf(temp, "%lf", val_arr[val_i - off]);
        } else if (equ[equ_i] == 's') {
          for (i = 1, sum = 0; i <= val_i - off; sum += val_arr[i], i++);
          sprintf(temp, "%lf", sum);
        }
        replace(equ + equ_i, equ+brac_i+1 , temp);
      }
    }
    val_arr[val_i] = cal(equ);
  }
  return val_arr[val_i-1];
}

double limrec_advan(const char *equ_argu, double *val, int val_i) {
    double y,z,d;
    recur_advan(equ_argu,val,val_i,Inf_r);
    y = val[Inf_r];
    z = val[Inf_r - 1];
    if (d < (double)1000*Appr) {
    return y;
   }
   else {
     cur_err=lim_diffuse;
     return -1;
   }
}

double recur(const char *equ_argu, double *val, int val_i, int n) {
    //val has starters beginning with index 0, val_i stands for the number of starters
    //n is the index of the wanted number (starting with index 1)
    int i;
    double *val_arr=malloc((n*2)*sizeof(double));
    for (i=0;i<val_i;i++) {
        val_arr[i+1]=val[i];
    }
    return recur_advan(equ_argu,val_arr,val_i+1,n);
}

double sumrec(const char *equ_argu, double *val, int val_i, int n) {
    //val has starters beginning with index 0, val_i stands for the number of starters
    //n is the index of the wanted number (starting with index 1)
    int i;
    double *val_arr=malloc((n*2)*sizeof(double));
    for (i=0;i<val_i;i++) {
        val_arr[i+1]=val[i];
    }
    recur_advan(equ_argu,val_arr,val_i+1,n);
    return sum_t(val_arr,n);
}

double limrec(const char *equ_argu, double *val, int val_i) {
    int i;
    double *val_arr=malloc((Inf_r+16)*sizeof(double));
    for (i=0;i<val_i;i++) {
        val_arr[i+1]=val[i];
    }
    return  limrec_advan(equ_argu,val_arr,val_i+1);
}

